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We model the evolution of the magnetic fields of neutron stars as consisting of a long 
term power-law decay modulated by short term small amplitude oscillations. Our model 
predictions on the timing noise v of neutron stars agree well with the observed statistical 
properties and correlations of normal radio pulsars. Fitting the model predictions to the 
observed data, we found that their initial parameter implies their initial surface magnetic 
dipole magnetic field strength Bo ~ 5x 10 14 G when to = 0.4 yr and that the oscillations 
have amplitude K ~ 10 -8 to 10 -5 and period T on the order of years. For individual 
pulsars our model can effectively reduce their timing residuals, thus offering the potential 
of more sensitive detections of gravitational waves with pulsar timing arrays. Finally our 
model can also re-produce their observed correlation and oscillations of i>, as well as the 
"slow glitch" phenomenon. 
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1. Introduction 

Many studies on the possible magnetic field decay of neutron stars have been done 
previously, based on the observed statistics of their periods (P) and period deriva- 
tives (P) . Some of these studies relied on population synthesis that includes compli- 
cated observational selection effects, unknown initial parameters of neutron stars, 
and not completely understood models of radio emission of pulsar j^ ^^^^^ l^^ 
consequently no firm conclusion c an be d rawn on if and how the magnetic fields 
of neutron stars decays (see R,e f s |S11QL11J f or rev iews). Alternatively some other 
studies used the spin-down or characteristics ages, r c = (P — Pq)/2P (Pq is the 
initial spin period of the given neutron star), or r c = P/2P if P ^> Pq, as indica- 
tors of the true age s of neutron stars and found evidence for their dipole magnetic 
11 911 ^1141 

field deca y 1 1 I However, as we have shown recently, their spin-down ages are 
normally significantly larger than the ages of the supernova remnants physically 
associated with them, which in principle should be the unbiased age indicators of 
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these neutron st ar JISIISHZ1 jt has been shown that this age mismatch can be un- 
derstood if the magnetic fields of these neutron stars decay significantly over their 
life times, since the magnetic fiel d decay alter s the spin-down rate of a neutron star 
significantly, so that T SN r/t c < ll l5 | 16 | 17 | 18 | 

Pulsars are generally very stable natural clocks with observed steady pulses. 
However significant timing irregularies, i.e., unp redicted times of arrivals of pulses, 
exist for all pulsars studied so far (see RefPS for an extensive reviews of many 
previous studies on timing irregularies of pulsars). The timing irregularies of the 
first type are "glitches", i.e., sudden increases in spin rate followed by a period 
of relaxation; it has been found that the timing irregularies of young pulsars with 
t c < 10 5 years are dominated by the recovery from previous glitch events^. In 
many cases the neutron star recovers to the spin rate prior to the glitch event and 
thus the glitch event can be removed from the data satisfactorily without causing 
significant residuals over model predictions. However in some cases, glitches can 
cause permanent changes to both P and P of the neutron star, which cannot be 
removed from the data satisfactorily. These changes can be modeled as permanent 
increase of the surface dipole magnetic field of the neutron star; consequently some 
of these neutron stars may grow their surface magnetic field strength gradually 
this way and eventually have surface ma gne tic field strength comparable to that of 
magnetars over a life time of 10 5-6 years^l These are the first studies that linked 
the timing irregularies of pulsars with the long term evolution of magnetic fields of 
neutron stars. 

The timing irregularies of the second type have time scales of years to decades 
and thus are normally referred to as timing noise. Hobbs et al. carried out so far the 
most extensive study of the long term timing irregularies of 366 pulsars^l. Besides 
ruling out some timing noise models in terms of observational imperfections, random 
walks, and planetary companions, some of their main conclusions are: (1) Timing 
noise is widespread in pulsars and is inversely correlated with t c ; (2) Significant 
periodicities are seen in the timing noise of a few pulsars, but quasi-periodic features 
are widely observed; (3) The structures seen in the timing noise vary with data span, 
i.e., more quasi-period features are seen for longer data span and the magnitude of 
\v\ for shorter data span is much larger than that caused by magnetic braking of 
the neutron star; and (4) The numbers of negative and positive v are almost equal 
in the sample. 

Proper understanding of these important observations not only can lead to better 
understanding of the internal structures of neutron stars, but can also be important 
for using some millisecond pulsars for gravitational wave detections. In particular 
it is of high interests if some of the long term timing noise can be modeled and 
removed reliably from data, thus improving the sensitivity of gravitational wave 
detections substantially. It is extremely desirable if a physical model of neutron 
stars can catch and even predict the main features and long term variations of the 
timing noise of pulsars, since such a model will ensure that no signal of gravitational 
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waves is removed when modeling the observed timing noise. Some of the observed 
periodicities have been suggested as due to free precessions of neutron stard^. Some 
strong quasi-periodicities have been identified recently as due to abrupt changes 
of the magnetospheric regulation of neutron stars, perhaps due to varied particle 

IOOI 

emissions^ 1 . It was suggested that such varied magnetospheric particle emissions 
is also responsible for the observed long term timing noised Alternatively it has 
been suggested that the propagation of Tkachenko waves in a neutron star may 
modulate its rotational mome nt of iner tia, and thus produce the observed periodic 
or quasi-periodic timing nois cF^^^ . Because none of the above processes may 
permanently change the properties of neutron stars, we can collectively represent all 
the above modulations as some sort of oscillations of the observed surface magnetic 
held strengths of neutron stars, as shown below. 

In this work we first model the evolution of the dipole magnetic field of a neutron 
star as a long term power-law decay modulated by one or several components of 
oscillations. We then calculate the spin-down history of a neutron star within the 
magnetic braking model with the prescribed dipole magnetic field. The timing noise 
is then calculated and compared with data. We show that all main observed fea- 
tures of neutron star's timing noise can be reproduced satisfactorily. In the end we 
compare our model with other recently developed models of neutron star's timing 
noise. 



2. Evolutionary model of the dipole magnetic field of neutron stars 

Attributing the observed Trn r / t c <C 1 to the long term decay of the magnetic field 
of a neutron sta r^^^^ , we found that the magnetic field decay should follow 
a power-law form, B oc t~ a (a w 0.5), in good agree ment w ith the field decay 
process dominated by the ambipolar diffusion mechanisnPlIZHl w ith constant core 
temperatures^^. We further found that the core temperatures of magnetars, nor- 
mal radio puls ars and millisecond pulsars are approximately 10 8 K, 10 7 K, and 10 5 
K, respectively^^! The neutron star free precession, inferred from the observed 
periodic timing nois is phenomenologically equivalent to periodic oscillation of 

its apparent surface magnetic field B = *J B^ sin 2 x, where x is the angle between 
the magnetic axis and spin axis of the neutron and varies slightly as the neutron 
star precesses periodically. Similarly the observed strong quasi-periodicities of tim- 
ing noise can also be phenomenologically modeled as modulations of the observed 
surface magnetic fields of these neutron stars, which are obtained from the observed 
P and P of neutron stars. We therefore model the evolution of the dipole magnetic 
field of a neutron star as a long term power-law decay modulated by one or several 
components of oscillations, 

i? = B (f)- a (l + E^ sin (^ + 2Tr ^))' (!) 

where t is the neutron star age, to and Bq are the starting time and initial surface 
dipole magnetic field strength of the neutron star, Ki, fa and Ti are the amplitude, 
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phase and period of the oscillating magnetic field of the i-th component, respectively. 

Assuming pure magnetic dipole radiation as the braking mechanism for a pulsar's 
spin down, we have, 

im = J -^V\ (2) 

and 

v = -AB 2 v 3 , (3) 

where A — y 6cil ' is assumed to be a constant, B is its dipole magnetic field at its 
magnetic pole, R is its radius, / is its moment of inertia. We then have, 

v = -iAB 2 v 2 v - 2ABBv 3 . (4) 

Simple calculations show that the value of the first term in the rig ht p art of the 
above equation is much less than the observed v as reported in RefP^H Therefore 
the observed timing noise characterized by v may be dominated by the second term 
containing B. From Eq. (p}, we have, 

B = fl(t)(__+£_^-*cos(& + 27r-)). (5) 

Assuming Ki -C 1 and ignoring the first term in Eq. Q, we have, 

^-^ ( T"5I^r cos(</,i+ ^ )} ' (6) 

where the first and second terms in the right part of the above equation are con- 
tributed by the long term magnetic field decay and short term magnetic field oscilla- 
tions, respectively. The age of the pulsar (not to be confused with its characteristic 
or spin-down age r c ) is given, 

For a young pulsar with t < 10 6 yr, it is very likely that the first term in Eq. (j6]) 
dominates and thus it is generally anticipated that v > 0. For example, for t ~ 

5 x 10 3 yr, a ~ 0.5, P ~ 1 s and v ~ 10~ 13 Hz s _1 , the first term in Eq. ([6]) gives 
v ~ 10 -24 Hz s _1 s _1 , consistent with observations. On the other hand, the second 
term in Eq. ([6]) dominates for an old pulsar, and thus both positive and negative 
v should be observed with almost equal possibilities. This is in general agreement 
with observations^! 



3. Statistical properties of the spin evolution of pulsars 

In estimating v from observations, is, i>, and v are obtained by fitting the phases of 
all pulses observed to the third order of its Taylor expansion over a period of time, 

*(t) = {<P + vt + ht 2 + ht 3 ). (8) 

2 b 
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We can therefore estimate v for a pulsar from 

v~-2i,(Z±f), (9) 

where / = 2ir max(Ki/Ti). As shown later, Ki <C 1 and Tj is on the order of years, 
therefore / <C 1/t for young pulsars. Substitute Eq. (JTJ) (with a — 0.5) into Eq. (|9|), 
our model of power-law magnetic field decay in neutron stars predicts, 

*~9o^± 1 (-, (10) 

'c r 'c 

where go = (3.2 x 10 19 ) 2 /(4^o-Bq) is the initial parameter of pulsar's magnetic field. 
In Fig. Q] (left) we show the observed correlation of v ~ P/t% for young pulsars with 
v > and r c < 2 x 10 6 yr. Fitting this correlation with that predicted by Eq. (fTU|) 
with / = 0, we get go = 83.74, which is considered to be a constant for all pulsars 
in the rest of this paper, unless specified otherwise. Assuming Bo = 5 x 10 14 G, 
we have to = 0.39 yr, suggesting that these pulsars may be fast rotating magnetars 
when t = to. In Fig. Q] (right) we show that the observed and model calculated v 
are very consistent with each other for these young pulsars, with / = in Eq. (jlOp ; 
this is clear evidence for the power-law decay of neutron stars' magnetic fields, at 
least for young pulsars. 




Log 10 v (s" 3 ) Log 10 v (s" 3 ) 

Fig. 1. Left: Correlation of v ~ P/t£ for young pulsars with v > and r c < 2 X 10 6 yr. The solid 
line is the best fit with Eq. (1 1 I t (/ = 0), resulting in go = 83.74. Right: Correlation between the 
observed and model calculated v for these young pulsars, with / = in Eq. dlOI . 

In Fig. [2] we show the observed correlation of v ~ r c for all pulsars with r c < 
10 9 yr (i.e. recycled millisecond pulsars are not considered); the left and right panels 
displays all pulsars with v > and v < 0, respectively. Eq. (fTUf means that v is 
determined by both P and r c , given go and /. Since the observed P of these pulsars 
is narrowly distributed around a median value of 0.6 s, therefore the observed main 
correlation should be with r c , rather than P, consistent with datcP^. In order to 
make straightforward comparison between data and that predicted by Eq. fjlO[) . we 
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Fig. 2. Correlation of v ~ t c for all pulsars with r c < 10 9 yr, i.e. recycled millisecond pulsars are 
not considered. Left: pulsars with v > 0. Right: pulsars with v < 0. The solid and dashed lines are 
the model predicted correlation of v ~ r c by taking the '+' and ' — ' signs in Eq. dlOI . respectively; 
different values of / are labeled in the figure. 

simply take P = 0.6 s in Eq. (|10[) to calculate the model predicted correlation of 
v ~ r c by choosing several values of /; the calculations with the '+' and ' — ' signs in 
Eq. (|10p are shown as solid and dashed lines, respectively. It should be noted that 
the calculations with the ' — ' sign in Eq. (|10p can also produce positive is, and thus 
dashed lines appear in both panels. For / = 10~ 13 and T = 3 yr, we have K ss 10~ 6 . 
In Fig. [3]we show the comparison between the observed and model predicted v. It 
is clear that with just two parameters (go is fixed and / ranges between 10 -15 and 
10~ n ), the general properties of the long term spin evolution of all pulsars can 
be described successfully, except for the recycled millisecond pulsars which we will 
address separately. This can be considered as a strong support to our model of the 
long term (power-law decay) and short term (small amplitude oscillation) magnetic 
field evolution of neutron stars. 




-30 -29 -28 -27 -26 -25 -24 -23 -22 -30 -29 -28 -27 -26 -25 -24 -23 



Log 10 |v| (s- 3 ) Logjvl (s 3 ) 

Fig. 3. Comparison between the observed and model predicted v with / = 3 X 10 13 or 10 12 for 
v > or < 0, respectively; here the observed value of P of each pulsar is used in Eq. ( 1 1 D . 
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4. Timing noise of individual pulsars 

Eq. (|10j) with fixed go and / cannot describe accurately the timing noise of individual 
pulsars for at least the following reasons: (i) go and / may be different for different 
pulsars; (ii) multiple oscillating components may exist with different phases, periods 
and amplitudes; (iii) the observed v for each pulsar is determined by fitting the 
observed pulse phases with Eq. ((8]) over a certain observation time span T s ; (iv) 
the observed v correlates with T s , which cannot be predicted by Eq. (|10p : (v) in 
using Eq. ([8]) to fit data, all higher order terms are absorbed in the four terms in 
Eq. (jH) ; and (vi) substantial residuals are often found in in using Eq. ((HJ to fit data. 
None of the last three issues can be addressed with Eq. (fT0|) . We therefore should 
calculate numerically the history of a pulsar's pulse phase within our model 
and compare it with observations. 
The pulse phase change is given by 

rf$ = {v + ut)dt, (11) 

where v and i> can be calculated numerically by combining Eqs. ([1]) and ©• The 
entire history of for a pulsar can thus be calculated by integrating Eq. ([11]) 
numerically, given the values of these parameters to, ol, (pi, Ki and Tj. The 
calculated can then be fitted with Eq. (|SJ) to determine v and residuals. 

4.1. Periodic and quasi-periodic residuals of timing noise 

For many pulsars the residuals in fitting data with Eq. (|SJ) show period or quasi- 
periodic behaviori^. From the FFT of the residuals of B1540— OfP^I, we can deter- 
mine the phase and period of the main peak in its power spectrum, which are taken 
as 4> and T in Eq. ([I} for a single oscillation component. The only remaining param- 
eter is K. With different values of K, we repeat the above procedure to calculate 
<!>(£), and fit it with Eq. ((SJ to determine residuals, which are then compared with 
data until the minimum \ 2 of the match is found. Fig. |4] (left) shows the comparison 
between the observed and model predicted residuals with K = 7.41 x 10 -10 ; the rms 
of the residuals is reduced from 11.7 ms to 7 ms by introducing just one parameter 
K, after subtracting the best fit residuals from the observed ones. 

When several quasi-periodic components are found in the residuals, we take an 
iterative approach to find the phase and period of each component: (1) find 4>, T 
and the amplitude of the strongest component in the FFT; (2) remove this peri- 
odic component from the residuals; and (3) repeat the above until no significant 
component is found in the FFT of the residual after all identified components are 
subtracted. We then still use the \ 2 optimization method to find A", of each compo- 
nent. In Fig. |4] (right) we show an example for B1826-17, the FFT of the residuals 
of which has two prominent components'^ the rms of the residuals is reduced from 
33 ms to 14 ms in this case. 
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Fig. 4. Residuals after fitting Eq. JHJ to the observed and calculated pulse phases for B1540— 06 
(left panel) and B1826-17 (right panel), respectively. The dashed and solid lines in the top panels 
are from the observed and calculated phases, respectively. The bottom panels are their differences. 



4.2. Oscillating v and slow glitches 

We then try to model the timing noise of B0329+54, which shows v correlates 
with T s and even switches between positive and negative values; the parameters we 
choose are given in the caption of Fig. [5] the left panel of which shows a similar 
correlation of v ~ T s in Fig. 12 of RefM This shows that a single pulsar can 
produce v over almost the entire range of v for all pulsars, depending on the chosen 
observation time span T s ; this may be the main source of the observed wide scatter 
in Fig. [2j In Fig. [5] (right), the calculated correlation of v ~ T s is shown with 
K = 10 -7 , exhibiting clear oscillations of v between positive and negative values. 
This explains naturally the almost equal numbers of positive and negative values of 
v reported^. 



Positive 
Negative 
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Fig. 5. Left: Calculated correlation of v ~ T B for B0329+54, with one oscillation component of 
T = 15.5 yr and K = 10 — 8 . Right: The same as the left panel, except that K = 10~ 7 . 



We then increase K to K = 2.5 x 10 5 and show in Fig. [6] the calculated 
correlations of v ~ T s and v ~ T s . We find that sometimes i> switches to positive 
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values, quite similar to the recently found slow 

glitches^MJJ 

the observed change 

of sign of v is also consistent with our model prediction. However the calculated 
v becomes negative most of the times. This means that the phenomenon of slow 
glitches may just be a manifestation of unusually large amplitude oscillations of the 
apparent magnetic field of a neutron star, consistent with the observed unusually 
large rms of residuals for the two pulsars observed with slow glitches, i.e., 1284.3 ms 
and 613.7 ms for PSR J1825-0935 (B1822-09) and PSR J1835-1106, respectively. 




Fig. 6. Calculated correlations of v ~ T B (left panel) and v ~ T s (right panel) with the following 
parameters: T = 15.5 yr and K = 2.5 X 10 -5 . 



5. Conclusions and discussions 

Our previous work has found evidence that the magnetic fields of neutron stars 
have long term decays in a power-law form, agreeing with the prediction of the 
ambipolar diffusion mechanism. Through examining the observed (apparent) short 
term variations of the magnetic fields of neutron stars, we propose that the (appar- 
ent) magnetic fields of neutron stars have oscillations. We thus model the evolution 
of a neutron star's magnetic field as consisting of a long term power-law decay 
modulated by short term small amplitude oscillations. Our model predicted v, a 
parameter widely used to characterize the timing noise of pulsars, agrees well with 
the observed statistical properties and correlations of all pulsars (except those recy- 
cled millisecond pulsars). This can be considered as a strong support to our model 
of magnetic field evolution of neutron stars. 

Fitting the model predictions to the observed data, we found that their initial 
parameter implies their initial surface dipole magnetic field strength Bq ~ 5 x 10 14 G 
when to = 0.4 yr and that the oscillations have amplitude K ~ 10 -8 to 10~ 5 
and period T on the order of years. We then calculated the timing properties of 
individual pulsars with our model. We modeled the observed timing residuals of 
two pulsars and demonstrated that our model can effectively reduce their timing 
residuals, thus offering the potential of more sensitive detections of gravitational 
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waves with pulsars. Depending on the combination K and T, we can re-produce the 
observed correlation and oscillations of v, as well as the "slow glitch" phenomenon. 

We did not study the timing noise properties of the recycled millisecond pulsars 
with our model, because their initial properties may be substantially different from 
the normal radio pulsars we studied in this work. 
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